Three-dimensional vortex configurations in a rotating Bose Einstein condensate 
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We consider a rotating Bose-Einstein condensate in a harmonic trap and investigate numerically the behavior 
of the wave function which solves the Gross Pitaevskii equation. Following recent experiments [6], we study in 
detail the line of a single quantized vortex, which has a U or S shape. We find that a single vortex can lie only 
in the x — z or y — z plane. S type vortices exist for all values of the angular velocity Q while U vortices exist 
for fi sufficiently large. We compute the energy of the various configurations with several vortices and study 
the three-dimensional structure of vortices. 
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I. INTRODUCTION 

Several experimental groups have produced vortices in 
Bose Einstein condensates (BEC) [1-6]. One type of experi- 
ments consists in imposing a laser beam on the magnetic trap 
holding the atoms to create a harmonic anisotropic rotating 
potential. For a prolate trap, it has been observed [2, 3, 6] that 
when a single vortex exists, the vortex line is not straight along 
the axis of rotation, but bending. Theoretical works [7, 8] es- 
tablish a simpler expression of the Gross Pitaevskii energy that 
only depends on the vortex lines. In [8], it is proved that bend- 
ing occurs for prolate condensates, but not for oblate ones. 

Minimization algorithms [9, 10] have been used to compute 
local minima of the Gross Pitaevskii energy and provide an 
evidence of the bending in the same setting as the experiment. 
Bending (or U) vortices are described in detail, and multiple 
vortex configurations are addressed in these studies. 

Recently, the ENS group [6] has further studied configu- 
rations with a single vortex line. They have observed planar 
bent vortices (U) but also different configurations (S). They 
study the length of the line, its deviation from the center and 
its angular momentum. 

In this paper, motivated by the recent experiments at the 
ENS [6], we numerically look for local minimizers of the 
Gross Pitaevskii energy and we want to understand the various 
vortex configurations observed in the experimental setting: U 
vortices but also S vortices. We compute solutions with up 
to 4 vortices and describe their three-dimensional structure. 
Different solution branches are followed and the evolution of 
the corresponding energy and angular momentum are shown. 
The framework of this study is the case of a prolate conden- 
sate where bending is an important phenomenon. 

We consider a pure BEC of N atoms confined in a harmonic 
trapping potential rotating along the z axis at angular velocity 
fl. The equilibrium of the system corresponds to local minima 
of the Gross-Pitaevskii energy in the rotating frame 
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where g 3 D = 4Trh 2 a/m and the wave function <\> is normal- 
ized to unity J v \(f>\ 2 — 1. Here, for any complex quanti- 
ties u and v and their complex conjugates u and v, (u, v) = 
(uv + uv)/2. 

For numerical applications, it is more convenient to rescale 
the variables as follows: r = x/R, u(r) = R 3 ' 2 (j)(x.), where 
R = d/ ^fe and 
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In this scaling the Thomas-Fermi limit of u is 

p TF (r) = Po - (x 2 + a 2 y 2 + (3 2 z 2 ) 



il = n/(euj x ). 
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Then, we use the dimensionless energy introduced in [7] 

E(u) =H{u) -ClL z (u), (2) 
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L z (u)=i^u(y^-x^), (4) 
defined in the domain T> — {^(r) > 0} . 



A. Numerical method 

In the present study we compute critical points of E(u) by 
solving the norm-preserving imaginary time propagation of 
the corresponding equation: 

du 1 1 

— - -Ait + i(Ct x r).Vu = -^u(p Te - \u\ 2 ) + fi e u, (5) 

with u = on dT> and /i e the Lagrange multiplier for the 
norm constraint j v \u\ 2 — 1. A hybrid 3 steps Runge-Kutta- 
Crank-Nicolson scheme [11] is used to advance the equation 
in time: 
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where TL contains the remaining non-linear terms. The corre- 
sponding constants for every step (I = 1,2,3) are : 

<zi = 8/15, a 2 = 5/12, a 3 = 3/4, 
6j = 0, b 2 = -17/60, 6 3 = -5/12, 
ci = 8/15, c 2 = 2/15, c 3 = 1/3. 

The resulting semi-implicit scheme is second order time accu- 
rate and allows reasonably large time steps, making it appro- 
priate for long time integration. The large sparse matrix linear 
systems resulting from the implicit terms are solved by an al- 
ternating direction implicit (ADI) factorization technique. 

For the spatial discretization we use finite differences on 
a Cartesian uniform mesh with periodic boundary conditions 
in all directions. To accurately resolve sharp gradients of 
the variable in presence of vortices, low numerical dissipa- 
tion and very accurate schemes are required for the spatial 
derivatives. A sixth-order compact finite difference scheme 
[12] with spectral-like resolution was chosen to this end. 
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The constants a v , (3 V control, respectively, the curvature and 
the height of the vortex. 

We first focus on single vortex configurations and describe 
later multi vortex configurations. 



II. SINGLE VORTEX LINES 

We have observed three different types of single vortex con- 
figurations as shown in figure 1 : planar U vortices, planar S 
vortices and non-planar S vortices. The U vortices are the 
bent vortices computed in [9, 10] and theoretically studied 
in [7, 8]. They are global minimizers of the energy. The S 
configurations were observed experimentally very recently [6] 
and are only local minimizers of the energy. 



B. Physical and numerical parameters 

The values of constants in (5) are set to e = 0.02, a 
1 .Of). (3 = 0.067, corresponding to the experiments of the 
ENS group [3, 10] (m = 1.445 • 10" 26 fc ff , a = 5.8 • 10- u m, 
N = 1.4 • 10 5 and uo x — 1094s -1 ). The angular frequency 
SI will be varied from to the maximum value of 0.9w x , for 
which no deformation of the condensate has to be taken into 
account. 

Equation (5) is propagated in imaginary time until the evo- 
lution of the energy (2) has a gradient in time smaller than 
10~ 6 . For the considered range of fi, the numerical do- 
main is fixed to an elongated box (x,y,z) 6 [—0.6,0.6] x 
[-0.6,0.6] x [-8.5,8.5]. A refined grid using 72 x 72 x 
510 nodes is employed, which is sufficient to achieve grid- 
independence for all considered numerical experiments. 

Different initial conditions are used in to trigger single or 
multiple vortex configurations and follow the corresponding 
branches as fi is varied. The simplest initial condition as- 
sumes a steady-state solution u(x, y, z) — ^ Ptf(x, y, z) and 
is useful to study vortex-free configurations and their degen- 
eracy into multiple vortex configurations when increasing the 
value of SI. Initial conditions with vortices are obtained by 
superimposing to the steady-state a simple ansatz for the vor- 
tex. For example, an initial condition with a centered straight 
vortex of radius e is obtained by imposing 



u(x,y,z)=*fp^-u e , (7) 
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where (r, ip) are the polar coordinates in the (x, y) plane. The 
3D shape of the vortex can be easily modified by shifting the 
center r of the vortex in successive (x, y) planes; for in- 
stance, to obtain a planar S shape vortex, the following func- 




FIG. 1: Single vortex configurations in BEC: (a) U vortex, (b) pla- 
nar S vortex, (c) non-planar S vortex. Iso-surfaces of lowest density 
within the condensate. 



A. U vortex 




FIG. 2: Single U vortex configurations for fl/ui^ = 0.42 (a), 0.58 
(b), 0.78 (c). 

The U vortex is a planar vortex formed of 2 parts: the 
central part is a line which stays on the z axis and the outer 
part reaches the boundary of the condensate perpendicularly. 
When SI increases, the central straight part gets longer (figure 
2) and the angular momentum (L z ) increases to 1 (figure 3). 

The U vortex is obtained by starting the simulation with an 
initial condition containing a straight vortex away from the z 
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axis. In fact, the U vortex lies either in the x — z or y — z 
plane. Starting with an initial condition which is not in one of 
these plane yields a final state in the y — z plane, which is the 
plane closest to the z axis. 

The shape of the the U vortex and its preferred location in 
the y — z plane can be analyzed using the approximate energy 
derived in [7, 8]: setting the vortex free solution to energy, 
then the energy of a vortex line 7 can be approximated by 



Ptf dl — C51 



dz, 



(8) 



where C is a constant which depends on the experimental pa- 
rameters and p TF is given by (1). If 7 is not in the x — z or 
y — z plane, then one can construct small perturbations of 7 
that preserve p TF and lower the energy. This implies that 7 
cannot be a critical point of the energy because the gradient 
is not zero. Of course, if the ellipticity of the cross section is 
small, the gradient is small, which may allow to observe these 
configurations. 

In order to understand the existence of the straight central 
part of the U vortex, one can also refer to the analysis of 
[8]: from equation (8) we can infer that a vortex line with 
a lower energy than the vortex free solution is obtained when 
the quantity p TF — C51p FF is negative, i.e. CTip TF > 1. Let 51 be 
such that CtipQ — 1; it corresponds to the 2d critical velocity 
for the existence of a vortex in the plane z — 0. For 51 close 
to 51, the inner region where CSlp TF > 1 is concentrated near 
the center of the condensate. In this region, the vortex line has 
to be straight (see [8]). This straight part is getting longer as 
51 increases since the region where CUp^ > 1 is getting big- 
ger. This region corresponds to Q > fl2d{z), where Q,2d{z) 
is the critical velocity for the existence of a vortex in the 2 
dimensional section where z is constant. In the outer region, 
the vortex reaches the boundary using the shortest path. 

Figure 3 shows the energy and angular momentum variation 
with 51 for the single vortex configurations. The U vortices 
exist only for 51 bigger than a critical value 51 c = 0.42^. It 
is interesting to note that at 51 c , the energy of the U vortex is 
bigger than the energy of the vortex free solution (we have set 
to zero the energy of the vortex free solution). A zoom in this 
region shows that Sl c is very close to the angular velocity 51 1 
for which the energy of the vortex free solution is equal to the 
energy of the U vortex. Figure 3 also shows that the angular 
momentum L z of the U vortex for 51 = 51 c does not go to 0. 
This suggests that in fact there could be another U solution 
for 51 > 51 c . Using an ansatz, another type of U solution is 
obtained in [10] which is a saddle point of the energy: it is 
away from the axis and has lower angular momentum. In [8], 
it is proved rigorously that for 51 small, there is no U as a 
critical point of the energy. 

For an initial condition with a straight vortex centered on 
the z axis, if 51 < 0.8w x , the straight vortex is unstable and 
the final configuration is a U, but if 51 > 0.8u> x , the straight 
vortex is stable. This is in agreement with the result of [8] 
where the local stability of the straight vortex for 51 larger is 
proved. 

For small 51, the U vortex disappears and a vortex-free con- 
figuration is obtained, while for 51 larger the U vortex degen- 
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FIG. 3: Energy (in units of Uui x ) and angular momentum per particle 
(in units of Ti) for the single vortex configurations. 



erates into a three-vortex configuration (described later). 



B. S vortex 

Motivated by the experiments of [6], we compute new criti- 
cal points of the energy, which are S configurations (see figure 
1). Several numerical experiments were performed, starting 
from different initial conditions containing an ansatz for the S 
vortex (see section IB). 

The planar S can be regarded as a U, with the half-part 
in the plane z < rotated with respect to the z axis by 180 
degrees (see figure 4). The non planar S are such that the 
projections of the branches on the x — y plane are orthogonal, 
i.e. the rotation of the branches is of 90 degrees. We could 
check that non planar S configurations with an angle between 
the branches different from 90 degrees do not exist. 




FIG. 4: Comparison between the single vortex configurations ob- 
tained for the same angular velocity fl/ui x = 0.44. Superposition of 
the U and S vortex (a) and the planar and non-planar S vortex (b). 

As already mentioned for the U vortex, stable planar S con- 
figurations lie either in the x — z or y — z plane. As for the 
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U, this can be explained using the limiting energy obtained in 
[8] and considering separately the upper or lower part of the 
S. As soon as the cross section is not a disc, if the upper or 
lower branch of the S configuration does not lie in the x — z 
or y — z plane, then the gradient of vortex line energy (8) can 
never be zero when 7 is varied. 




FIG. 5: Single S vortex configuration for Q /lj x = 0.38 (a), 0.44 (b), 
0.48 (c). 
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FIG. 6: Single vortex configuration. 



The S vortices exist for all values of while the U only 
exist for ft > tt c . When il decreases, the extension of the 
S along the z axis goes downwards as shown in figure 5, the 
angular momentum decreases to (figure 3) and the vortex 
tends to the horizontal axis. Note that a vortex along the hor- 
izontal axis has L z = 0, but a positive energy. On the other 
side, when increases, the S gets straighter and it tends to 
the vertical axis. 

The global minimum of the energy is never an S. But the 
difference in energy (and angular momentum) between U and 
S vortices is very small, as illustrated in figure 3 because an 
S vortex is almost like a U with a half-part rotated by 180 
degrees. 



C. Minimizer with fixed L 

As pointed out in [6], the minimization problem which is 
related to the experiments, is rather to minimize H (see (3) 
while fixing L z , rather than minimizing E = H — QL Z . This 
has been studied in the 2 dimensionnal setting in[13]. One can 
notice that if a given configuration with H = h and L z — I 
minimizes E = H — ttL z for some f2, then h minimizes H 
under the constraint that L z = I: indeed if W = H(u) with 
L z (u) = I, then H' — 0,1 > h — Ql, since (h, I) minimizes E, 
and this implies that H' > h. Moreover fl is the slope to the 
curve H(L Z ) at the point (h, I) and the property of minimizing 
E that is for all ti, V, 

ti - 9,1' > h - fll 

implies that the curve H(L Z ) lies above its tangent at this 
point. 

We have plotted H as a function of L z . We can check that 
the curve is convex, and above its tangent, which is consistent 
with the fact that we have computed minimizers of the energy. 

We know that the U solution exists for O > O c and has 
L z > 0.4. For L z < 0.4, we expect that the process of min- 
imizing H with fixed L z would produce U vortices and the 
curve H(L Z ) should be concave in this region. In [8], we 



have proved that for L z close to 0, H > CL 
first indication to the concavity of the curve 
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IH. MULTIPLE VORTICES 

Multiple vortex configurations are obtained based upon dif- 
ferent numerical strategies. The first one is to start the compu- 
tation from a vortex-free steady state and to abruptly increase 
O to a very hight value; multiple vortices are thus obtained. 
The second strategy is to generate an initial condition with 
vortices as described in section I B (the advantage being the 
control of the shape and initial arrangement of the vortices). 

Both techniques are used to follow solution branches with 
two, three or four vortices in the condensate. Figures 7 and 
8 display energy and angular momentum vs fl for all studied 
configurations. 
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FIG. 7: Energy (in units of huj x ) for all studied configurations. 
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FIG. 8: Angular momentum L z (in units of 7i) for all studied config- 
urations. 



A. 3 vortices 

When £1 is increased, the single U vortex solution switches 
to a 3 vortex configuration (£1 = 0.9u) x ), As shown in figure 
9a, the configuration is invariant under rotation in a central 
plane near z — but not near the edges. For large £1, three- 
dimensional views show (figure 9 a,b) that there are 2 vor- 
tices of similar size and a longer one which is bending near 
the boundary. For fi = 0.8uj x , all vortices display contorted 
shapes (figure 9c), very similar to those reported in [9]. Let 
us point out that the angular momentum of all these 3 vortex 
configurations is lower than 3 (see figure 8). 




FIG. 9: Three- vortex configuration for Q,/u x = 0.9 (a), 0.72 (b), 
0.68 (c). Lower pictures show iso-contours of |u| in the central z = 
cut plane. 

When we put as initial condition a configuration with 3 
identical U vortices at 120°, in the final state, one of them gets 
a little longer (figure 10a) and the symmetry is lost. This con- 
figuration has almost the same energy and angular momentum 
as the configuration displayed in figure 9b. In exchange, the 
initial condition with three straight vortices on the x-axis has 
its symmetry preserved (figure 10 b), but with a higher energy 






FIG. 10: Three-vortex configuration obtained for the same Q/u) x = 
0.72, from different initial conditions: 3 identical U vortices at 120° 
(a) and 3 straight vortices in a row on the x-axis (b). Lower pictures 
show iso-contours of |u| in the central z — cut plane. 

When further decreasing f2, the 3 vortex branch switches to 
a 2-vortex displaying irregular shapes (figure 1 1). 





FIG. 1 1 : Two vortices obtained from the 3-vortex configuration when 
the value of Q,/lo x is decreased to 0.64 (a) and 0.6 (b). 



B. 2 vortices 

The two-vortex branch presented in this section was ob- 
tained by starting from a vortex-free solution and suddenly 
increasing £1 to a value of 0.8^. The configuration is planar 
and symmetric, like twice a single U vortex, but away from 
the axis (there is a repulsion between the lines). 

When £1 increases, the lines are almost straight and get 
closer to each other. This is in agreement with the fact that 
when £1 gets large, the straight vortex is a local minimizer of 
the energy. Hence the bending is no longer the important phe- 
nomenon. 

We recall that decreasing £1 from a configuration with 3 
vortices, we obtained 2 vortices which are not symmetric, one 
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FIG. 12: Configuration with two symmetric vortices for Q./u x 
0.48 (a), 0.6, (b) 0.8 (c). 



FIG. 13: Four- vortex configurations for (a) Q/lj x = 0.86 - obtained 
from an initial condition without vortices and (b) Q/cu x — 0.72 - 
obtained from an initial condition with four symmetrical vortices. 



We have to point out that with the initial condition of 4 
identical vortices, the symmetry is preserved as displayed in 
figure 13b, which was not the case for 3 vortices. 



being longer than the other (figure 11). This configuration has 
slightly bigger energy than the 2 symmetric vortices. 



C. 4 vortices 



Starting from an initial condition without vortices and in- 
creasing to 0.86cj x , we have obtained stable configurations 
with 4 curved vortices (figure 13 a). When decreasing f2, 
this configuration rapidly degenerates into a three-vortex state. 
For lower we could obtain stable configurations with four 
symmetric vortices (figure 13 b), but with bigger energy. The 
location of the vortices in the plane z — is the same. 
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